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^1.0  INTRODUCTION 

^“Linear  elastic  fracture  mechanics  has  gained  a  substantial 
acceptance  in  industry  and  has  become  one  of  the  most  important 
design  considerations.  It  Is  now  well  recognized  that^nany 
engineering  structures  such  as  airplanes,  turbines,  piping  (pres¬ 
sure  vessels),  bridges,  etc.  contain  pre-existing  flaws.  As  a 
result  of  even  rather  moderate  service  loads,  crack  propagation 
resulting  from  these  flaws  can  have  a  dramatic  effect  on  the 
service  life  of  the  component.  To  account  for  this  reduction  in 
service  life,  fracture  mechanics  analysis  in  conjunction  with  a 
fracture  control  plan  is  generally  implemented. 

The  basic  elements  of  a  fracture  control  plan  have  been 
described  by  Rolfe  and  Bansom  (Ref.  1)  as  follows: 

•  **'■  w 

L  l  '  1.  Identif icatlon  of  the  factors  that  may  contribute  to  the 

structure.  Description  of  service  conditions  and 
loadings . 

2.  Establishment  of  the  relative  contribution  of  each  of 
these  factors  to  a  possible  fracture  In  a  member. 

j.  Determination  of  the  relative  efficiency  and  trade-offs 
of  various  design  methods  to  minimize  the  possibility  of 
failure . 

4.  Recommendation  of  specific  design  considerations  to 
ensure  the  safety  and  reliability  of  the  structure 
against  fracture. 

The  life  of  the  structural  component  is  generally  determined 
by  the  time  necessary  to  initiate  a  crack  and  to  propagate  the 
crack  from  a  sub-critical  to  critical  size.  Two  parameters  are 
required  for  successful  determination:  the  fracture  behavior  or 
the  material  and  the  state  of  stress  and  strain  around  the  crack 
tip.  The  three  measures  of  the  severity  of  stresses  and  strains 
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around  the  tip  of  cracks,  typically  employed  in  linear  fracture 
mechanics,  are  the  elastic  stress  intensity  factors  Kj,  Kji*  and 
Kjjj  for  the  opening,  the  inplane  shear  and  the  anti-plane  shear 
modes,  respectively. 

Numerous  methods  are  now  available  for  determining  stress 
intensity  factors.  Tada,  Paris,  and  Irwin  (Ref.  2)  present  a 
variety  of  methods  to  predict  these  factors,  including:  boundary 
collocation,  successive  boundary  stress  correction,  and  finite 
element  methods.  The  most  powerful  method  of  the  three  i3  the 
finite  element  method.  A  large  number  of  papers  have  been 
written  on  this  subject  alone.  Most  of  these  are  restricted  to 
two  dimensional  methods  and  linear  elastic  materials. 

Many  of  the  papers  written  on  finite  elements  used  for  frac¬ 
ture  studies  are  classified  as  either  hybrid  or  singular  element 
formulations.  Many  of  the  elements  developed  suffered  from  either 
lack  of  accuracy,  generality,  or  consistency.  Barsoum  (Ref.  3) 
points  out  shortcomings  of  several  different  elements.  These 
shortcomings  include  inability  to  model  rigid  body  or  constant 
strain  modes,  inability  to  include  thermal  or  body  force  effects, 
and  lack  of  compatibility  with  other  elements. 

The  discussion  presented  herein  attempts  to  illustrate  these 
shortcomings  and  suggests  alternative  two-dimensional  crack  ele¬ 
ment  formulations  that  are  less  restrictive.  This  report  presents 
the  theory,  implementation,  instructions,  and  sample  problems 
which  will  allow  COSMIG/NASTRAN  users  to  utilize  the  developed 
crack  element.  The  contents  of  this  report  are  as  follows: 

Section  2  presents  the  theoretical  development  of  the  crack  ele¬ 
ment;  Section  3  describes  how  the  crack  element  was  implemented 
into  COSMIC/NASTRAN;  Section  ^4  presents  various  numerical  results 
obtained  using  the  crack  element  and  Section  5  presents  detailed 
user  Instructions  and  sample  problems.  A  summary  of  the  work 
performed  and  conclusions  are  presented  in  Section  6. 
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2.0  THEORETICAL  DEVELOPMENT 


An  early  concept  for  predicting  crack  propagation  comes  from 
the  work  of  Griffith  (Ref.  4)  on  the  fracture  of  glass.  The 
basic  idea  of  his  theory  being  that  the  surface  of  a  solid  pos¬ 
sesses  surface  tension,  similar  to  liquids;  thus,  when  a  crack  in 
a  solid  propagates,  the  increase  in  externally  added  or  internal¬ 
ly  released  energy  is  balanced  by  the  increase  in  surface  tension 
energy.  If  in  an  elastic  solid,  V  and  U  represent,  respectively, 
the  work  of  the  externally  applied  forces  and  the  strain  energy, 
and  if  the  specific  surface  tension  energy  is  denoted  by  y»  then 
Griffith's  energy  balance  criterion,  as  shown  by  Reference  5,  may 
be  expressed  as  follows: 


(1) 


It  is  pointed  out  that  the  foregoing  energy-balance  relation 
is  only  a  necessary  condition  for  crack  growth,  with  the  quantity 
on  the  left-hand  side  representing  the  energy  available  for  frac¬ 
ture,  and  the  quantity  on  the  right  hand  side,  the  resistance  of 
the  lid  to  fracture  propagation.  Prom  this  physical  meaning  of 
the  terms  involved  in  the  energy-balance  equation,  it  also  follows 
that  the  stability  of  quasi-fracture  propagation  may  be  determined 
from : 


>  0  :  unstable  crack  growth 

3X  “  U)  -  y]  =  0  :  neutral  equilibrium  (2) 

<  0  :  stable  crack  growth 

Based  on  the  asymptotic  solutions  to  crack  problems  presented 
by  Westergaard  and  Sneddon  (Refs.  6  and  7),  Irwin  (Ref.  8)  gave 
the  name  of  "stress  Intensity  factor"  to  the  coefficient  which 
appears  in  the  asymptotic  expression  for  stress.  Irwin  noted 
that  the  energy  available  for  fracture  per  unit  crack  extension 
may  be  directly  related  to  that  coefficient  by: 


Where  the  quantity  G  (after  Griffith),  Introduced  by  Irwin,  Is 
known  as  the  "strain  energy  release  rate."  Subsequently,  Irwin 
showed  that  the  stress  and  displacement  fields  around  a  crack  tip 
In  a  linearly  elastic  solid  under  the  most  general  loading  condi¬ 
tions  may  be  expressed  in  terms  of  three  stress  intensity  factors 
Ki,  K2,  and  {(3.  associated,  respectively,  with  the  opening,  in¬ 
plane  shear  and  anti-plane  shear  modes  of  deformation. 

Generalizing  Irwin's  findings  for  linear  materials,  the 
asymptotic  expressions  for  stresses  and  displacements  near  the 
tip  of  a  crack  have  the  following  form: 


k  { C )  k  (t) 

°ij  ‘  -<r-  rij(0-ck>  *  —  ‘V9-‘V 

r  r 


ui  "  ET  k^tjr1"®  pJ(0,Ck)  +  k2(t)r1  “  P^(e,ck) 


for:  0  <  a  <  1  ;  1,J  =  x,y 


k  (t) 

aiz  “  flz(9’Ck)  ;  0  <  B  <  1  ;  1  -  x,y 

r 


uz  "  Tpr  k 3 ( t ) r 1  3  P3<  0,Ck) 


where  E*  and  p*  are  normalizing  material  moduli;  the  Ck  are 
dimensionless  material  constants;  and  f^  and  are  known, 
bounded  functions.  Also,  for  clarity,  the  following  notation  has 
been  used: 


In  the  previous  expressions,  the  powers  a  and  0  of  the 
r Lngularltles  differ  from  1/2  only  in  nonhomogeneous  and  in 
certain  homogeneous  but  anisotropic  materials. 

Quite  naturally,  in  dynamic  problems,  the  stress  intensity 

factors,  are  functions  of  time  and  are  defined,  following 
Irwin,  as  the  coefficients  of  the  singular  terms  in  the  expres¬ 
sions  for  stresses;  thus,  for  Mode  I  type  of  loading,  for  example 


def 

(t)  = 


lim  r 


f  (0,C.  )  r +0 

yy  *  k 


(6) 


Practical  use  of  the  previous  concepts  is  possible  when  the 
resistance  to  fracture  of  the  material  is  known.  Hence,  if  in 
Mode  I  fracture,  G  is  used  to  characterize  the  material,  the 
necessary  condition  for  fracture  becomes: 


G  i  =  G. 

1  lc 


(7) 


where  ,  the  critical  resistance  parameter,  is  known  as  the 
"critical  strain  energy  release  rate"  in  Mode  I.  Similarly,  if 


the  corresponding  critical  value  of  is  used  to  represent  the 


material's  resistance  to  fracture,  the  fracture  criterion  becomes 


K  i  =  K. 

1  lc 


(8) 


in  any  event,  since  G,  and  K,  . 

lc  lc* 


being  material  parameters,  are 
constant,  the  stability  of  crack  growth  would  be  determined  from 
dG1/da  and  dK-^/da,  respectively. 


More  generally,  under  three-dimensional  loading  conditions 
all  three  modes  of  deformation  mentioned  above  are  present  and 


since  G  Is  a  scalar  quantity,  the  total  energy  available  for 
fracture  Is  given  by: 


G  =  +  Gg  +  G3  (9) 

with  the  new  Incremental  crack  surface,  dA,  lying  In  the  plane 
that  corresponds  to  the  maximum  available  energy,  G,  provided  the 
material  Is  Isotropic  with  regard  to  fracture  resistance. 

In  addition,  fracture  criteria  in  terms  of  stress  intensity 
factors  usually  adopt  the  form  of  Interaction  envelopes  as  shown 
by  References  9  and  10,  either 


Ki  2 

0  ♦  ( 
lc 


K0  2  K,  2 

*  (it-) 

2c  *3c 


(10) 


or 


a  1 1  Kl  +  2al2KlK2  *  a22KI  *  a33K3  *  1  (U) 

as  suggested  by  Erdogan  and  Slh  (Ref.  11)  and  by  Slh  (Ref.  12) 
based  on  energy  arguments. 

Based  on  the  foregoing  presentation,  it  is  only  natural  that 
two  basic  approaches  have  been  followed  to  ascertain  the  fracture 
behavior  of  linear  elastic  solids  containing  cracks: 

1.  determination  of  stress  Intensity  factors 

2.  determination  of  strain  energy  release  rate 

Both  methods  are  essentially  equivalent. 

For  the  present  case,  the  determination  of  the  stress  inten¬ 
sity  factors  for  a  two  dimensional  system  using  finite  element 
formulation  will  be  considered.  Other  authors  have  suggested 
some  of  the  shortcomings  of  earlier  finite  elements  developed  for 
fracture  mechanics  studies.  The  elements  developed  by  Barsoum 
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(Ref.  3)  and  Henshell  and  Shaw  (Ref.  13)  rectified  many  of  the 
problems  described  above;  however,  these  elements  were  limited  to 
displacement  of  the  form  r^/2.  Consequently,  they  could  only 
model  strain  singularities  of  the  form  r“l/2.  Recently,  Stern 
(Ref.  1*4)  and  more  recently,  Hughes  and  Akin  (Ref.  15)  introduced 
families  of  consistent,  conforming  elements  which  allow  displace¬ 
ments  of  the  form  rT.  While  the  Stern  element  appears  to  have 


the  restriction  that  0  <  y  <  1»  the  element  of  Hughes  and  Akin  is 
valid  for  all  y  >  0.  The  element  described  herein  is  based  upon 
shape  functions  suggested  by  Hughes  and  Akin. 

The  element  presented  here  possesses  the  required  rigid  body 
and  constant  strain  modes.  It  properly  models  thermal,  body 
force,  and  pressure  loading  conditions.  Additionally,  it  is 
compatible  with  standard  linear  or  quadratic  isoparametric  ele¬ 
ments  and  can  be  used  as  a  nonsingular  element  with  a  variable 
number  of  nodes. 

The  following  sections  present  the  element  formulation. 

This  includes  the  assumed  shape  functions,  procedures  for  cal- 
culatirg  the  stiffness  and  mass  matrices,  and  equivalent  thermal 
loads.  Also,  the  equations  used  to  evaluate  stresses  and  stress 
Intensity  factors  are  presented. 

2.1  ELEMENT  FORMULATION 

The  following  derivation  follows  Hughes  and  Akin.  Referring 
to  Figure  1,  the  standard  bilinear  shape  functions  are  used  for 
grids  1  through  4: 

N1(r,s)  =  (l-r)(l-s) 

N2(r,s)  =  r(l-s) 

N3(r,s)  =  rs 

Nli  ( r , s  )  *  (l-r)s  (12) 


% 
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The  shape  functions  for  grids  5-8  are  chosen  as: 


N5(r,s)  =  (l-s)P(r,Y) 
NgCr.s)  =  rP(s,y) 

N7 ( r, s )  =  sP(r, y) 

Nq ( r , s  )  =  (l-r)P(s , y) 

where 

p(x.y)  -  2(X  - 

1  -  2 (1/2  )  Y 


(13) 


(14) 


It  can  be  easily  shown  that  the  shape  functions  for  grids 
5-8  reduce  to  the  standard  quadratic  serendipity  element  when 
Y  of  Equation  (14)  is  set  equal  to  2.  It  can  also  be  observed 
that  the  shape  function  for  grids  5-8  satisfies  the  interpola¬ 
tion  property  at  all  nodes  of  the  element.  That  is: 


W  ‘  5ij  and  W  =  hj 

where  rj  and  Sj  are  values  of  r  and  s  at  grid  j,  and  6^  is  the 
Kronecker  delta.  However,,  the  shape  functions  associated  with 
grids  1  -  4  do  not  satisfy  the  interpolation  property  at  grids 
5-8.  Following  the  standard  technique  (Ref.  15)  the  shape 
functions  for  grids  1-4  are  modified  as  follows: 


Nx  «-  N][(r,s) 
N2  +  N2(r,s) 

N3  *  N3(r»s) 
N4  ♦  N4(r,s) 


[Ng(r,s)  +  N^(r,s)  ]/2 
[Np.  (  r,  s )  +  Ng  ( r,  s  )  ]/2 
[Ng(r,s)  +  ( r,  s  )  ]/2 

[N  (r,s)  +  Ng(r,s) ]/2 


(15) 


where  the  <-  reads:  "is  replaced  by". 

It  can  be  shown  that  the  shape  functions  for  all  eight  grids 
satisfy  the  required  interpolation  property.  Additionally,  the 
shape  functions  are  capable  of  exactly  representing  the  monomials 
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1,  r,  s,  rY,  rs,  s2,  rYs,  and  s2r.  The  presence  of  1,  r,  and  3 
ensure  representation  of  rigid  body  and  constant  strain  modes. 

The  presence  of  rY  allows  exact  representation  of  displacements 
of  the  form  rY.  Note  that  this  will  result  in  a  line  singularity 
of  the  form  r^-^  upon  differentiation. 

In  order  to  represent  point  sigularities ,  the  quadrilateral 
form  is  degenerated  into  a  triangle.  This  is  done  by  coalescing 
grids  4,  8,  and  1  as  can  be  done  for  standard  isoparametric  ele¬ 
ments  (Ref.  16)  and  as  is  shown  schematically  in  Figure  2.  Thus, 
for  a  point  singularity,  the  shape  function  associated  with 
grid  1  is  replaced  with: 

N-^r.s)  ♦  N1(r,s)  +  N^r.s)  +  Ng(r,s)  (16) 

In  summary,  for  the  6-grld  triangle,  the  shape  function 
associated  with  grid  1  is  given  by  Equation  (16),  the  shape  func¬ 
tions  associated  with  grids  2  and  3  are  given  by  N2  and  N3  of 
Equation  (15),  and  the  shape  functions  associated  with  grids  5 
through  7  are  given  by  N5  through  N7  of  Equation  (13)* 

When  the  6-grld  triangle  is  used  and  y  of  Equation  (14)  is 
set  appropriately,  then  a  singular  element  or  crack  element  is 
developed.  If  y  is  set  equal  to  2,  then  the  standard  serendipity 
element  (Ref.  17)  is  obtained.  Additionally,  if  y  =  2,  any  of 
the  mid-side  grids  may  be  omitted.  This  element  can  then  be  used 
as  a  transition  element  to  change  from  a  mesh  of  quadratic  ele¬ 
ments  to  one  of  linear  elements.  If  the  singular  triangular 
element  is  used  (Figure  2),  then  the  only  grid  which  may  be  elim¬ 
inated  from  Figure  2  is  grid  6.  User  instructions  for  eliminating 
the  mid-3ide  grids  will  be  presented  in  Section  5« 

2.2  STIFFNESS  MATRIX,  MASS  MATRIX  AND  THERMAL  LOAD  VECTOR 

Given  the  shape  functions  of  the  previous  section,  calcula¬ 
tion  of  the  element's  stiffness  matrix,  mass  matrix,  and  thermal 
load  vector  follows  the  standard  procedure  as  described  in 
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Reference  17.  These  quantities  are  given  in  terms  of  element 
coordinates  as 


K6  =  /  BT  D  B  dV 

•  a,  /v 


MS  =  /  p  NT  N  dV 

ve 


P®  =  /  BT  D  a  AT  dV 


where 


B  =  L  N  ,  H  ■  [Nj  I,  M2  I, 


,  k  “ 


fx  0 


0  ay 

3  3 

ly  lx 


I  is  a  2  by  2  identity  matrix  and  the  definition  of  D  will  be 
given  in  Section  2.3.  See  Reference  17  for  more  details.  The 
integrations  are  performed  using  Gaussian  quadrature.  That  is, 
the  integrals  are  approximated  as: 


/  f(x)  dx  *  1  W.  f (a  ) 


(18) 


Due  to  the  formulation,  it  can  be  shown  that  along  the 
s  direction,  the  integration  order  needs  to  be,  at  most,  A  to 
exactly  integrate  the  element.  For  singular  forms  of  the  element, 
the  Integration  order  recommended  along  the  s  direction  is  4  and 
along  the  r  direction  is  5.  If  the  element  is  used  in  a  non¬ 
singular  form,  then  2  by  2  integration  is  recommended  for  undis¬ 
torted  or  slightly  distorted  elements,  and  3  by  3  integration  is 
recommended  for  distorted  elements  (Ref.  16). 


1  2 


2.3  STRESS  AND  STRESS  INTENSITY  FACTOR  CALCULATIONS 

For  the  present  element,  stresses  are  calculated  and  reported 
at  the  natural  coordinate  centroids.  These  correspond  to  the 
locations  of  s  =  r  =  1/2  in  Figures  1  and  2.  The  stresses  are 
calculated  using  the  equations 

o  =  De  =  DBu6 

where 


ue  are  the  element's  grid  displacements,  and  B  was  defined  In 
Section  2.2.  This  D  matrix  Is  for  plane  strain.  For  plane 
stress,  X  is  replaced  by  2\y/(A+2y). 

In  addition  to  the  stresses,  the  element  coordinates  of 
these  stress  locations  are  also  reported.  These  x  and  y  locations 
are  measured  in  element  coordinates  as  depicted  in  Figure  3»  The 
x  and  y  locations  are  given  by: 

x  =  I  N  x® 

1  X  1 


y  58  l  N 

i 


i  y? 


where  x|  and  y®  are  the  coordinates  of  the  element's  grid  points 
measured  in  element  coordinates  and  the  summation  is  carried  out 


over  the  number  of  grid  points.  The  shape  functions,  N^»  are 
evaluated  at  s  =  r  =  1/2. 
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The  stress  intensity  factors  calculated  are  based  upon  the 
stress  or  displacement  fields  of  the  crack  element.  When  the 
calculations  are  based  on  stresses,  the  resulting  equations  are: 

1/2 

Kt  =  lim  ( 2 irr)  '  o  ( 0=0) 
r+0  y 

Ktt  =  lim  (2*r)1/2  t  ( 0=0 )  (19) 

1  r+0  xy 

where  the  nomenclature  is  shown  in  Figure  4.  If  the  stress  in¬ 
tensity  factors  are  based  on  displacements,  the  equations  used 
are : 

Ki  =  yrr^r  ^)1/2  uy  (9-’) 

kii  ’  ^ 2  ux  <9"")  <20) 

Equation  (20)  is  for  plane  strain.  For  plane  stress,  v  is 
replaced  by  v/(l+v). 

The  limits  of  Equations  (19)  and  (20)  are  determined  by 
evaluating  the  expression  on  the  right  hand  3ide  of  the  limit 
3ign  at  the  Gauss  integration  points  along  the  ray  nearest  the 
grid  1-2  edge  depicted  in  Figure  3-  The  expressions  on  the  right 
hand  side  of  the  limit  signs  are  then  extrapolated  to  r  =  0  using 
Lagrangian  Interpolation. 

If  the  stress  intensities  are  based  upon  displacements 
(Equation  (20)),  then  the  crack  element  configuration  must  be  as 
shown  for  Element  1  in  Figure  5a.  If  the  stress  intensities  are 
based  on  stresses  (Equation  (19)),  the  crack  element  configura¬ 
tion  must  be  as  shown  for  Element  *4  in  Figure  5b.  Further 
instructions  regarding  calculations  of  stress  intensities  are 
presented  in  Section  5. 


(a)  Stress  intensity  factors  for  E 
be  based  on  displacements  (IKI 


3.0  IMPLEMENTATION  INTO  NASTRAN 

Implementation  of  the  crack  element  into  NASTRAN  was  per¬ 
formed  via  the  dummy  element  CDUM1.  This  procedure  is  discussed 
more  completely  in  Reference  18.  The  present  element  was  modeled 
after  the  CIS2D8  element  routines,  due  to  their  similarity. 

The  first  step  was  to  create  a  subroutine  KDUM1  which  gen¬ 
erates  the  stiffness  and  mass  matrices.  The  mass  matrix  may  be 
either  consistent  or  lumped.  When  the  mass  matrix  is  used  for 
calculation  of  gravity  loads,  the  consistent  mass  matrix  should 
be  specified.  This  subroutine  is  eventually  linked  to  NASTRAN 
LINK  8. 

For  computation  of  thermal  loads,  the  subroutine  EDTL  must 
be  modified  to  make  a  call  to  SSGETD  before  calling  the  routine 
DUM1.  The  dummy  coding  in  routine  DUM1  is  then  modified  to 
calculate  the  thermal  load  vector  based  on  the  average  element 
temperature.  After  EDTL  and  DUM1  have  been  modified,  they  must 
be  linked  to  NASTRAN  LINK  5. 

Finally,  the  dummy  coding  for  the  SDUM11  and  SDUM12  routines 
must  be  modified  so  that  It  performs  the  required  operations. 
SDUM11  performs  the  preliminary  geometry  calculations  and  creates 
various  data  arrays  to  be  used  in  Phase  2  stress  recovery.  SDUM12 
then  uses  these  data  arrays,  grid  point  displacements,  and  tem¬ 
peratures  to  compute  stresses  and  stress  intensity  factors  and 
writes  them  to  the  output  file.  After  the  SDUM11  and  SDUM12 
coding  has  been  modified,  it  is  linked  to  NASTRAN  LINK  13. 

An  overview  of  the  various  routines  and  their  functions  is 
presented  in  Figure  6.  These  subroutines  must  be  compiled  and 
linked  to  the  NASTRAN  executable  code.  Instructions  for  compiling 
and  generating  the  new  NASTRAN  executable  code  on  a  VAX  computer 
are  presented  in  the  appendix. 


4.0  NUMERICAL  RESULTS  AND  VERIFICATION  PROCEDURE 

To  assess  the  accuracy  of  the  present  element,  four  different 
crack  geometries/loading  conditions  with  known  solutions  were 
analyzed.  Figure  7  shows  the  different  geometries  analyzed. 

Figure  8  presents  four  different  mesh  sizes  which  were  used  to 
analyze  the  first  three  crack  geometries.  Figure  9  shows  the 
boundary  conditions  used.  For  the  edge  crack  with  a  point  load. 
Figure  9  is  modified  so  that  the  load  Is  applied  at  the  edge  of 
the  crack.  Table  1  presents  the  errors  associated  with  both  the 
crack  opening  displacement  (COD)  and  the  mode  I  stress  intensity 
factor  Kj .  As  can  be  seen,  the  COD  is  less  sensitive  to  the  mesh 
size,  while  the  Ki  values  appear  to  be  converging  to  their  exact 
solutions.  However,  the  edge  crack  with  a  point  load  solution 
appears  to  overshoot  the  exact  solution  by  about  5%*  It  should 
be  mentioned  that  the  "exact”  solution  for  the  edge  crack  speci¬ 
men  with  a  point  load  is  considered  to  be  accurate  to  within  2%. 
The  other  solutions  were  considered  to  have  accuracies  better 
than  1%.  These  solutions  were  obtained  from  Reference  2. 

Figure  10  presents  a  model  of  a  central  crack  In  a  finite 
plate.  To  ascertain  the  accuracy  of  the  element's  mode  II  stress 
intensity  factor,  the  model  of  Figure  10  was  used.  The 

results  for  both  Kj  and  Kji  are  presented  In  Table  2.  As  can  be 
seen,  the  Kjj  i3  within  about  4%  of  the  exact  solution. 

In  addition  to  the  test  cases  previously  described,  various 
simple  patch  tests  were  performed  to  ensure  proper  coding  of  the 
element  routines.  These  tests  include  rigid  body  motion  tests, 
free  thermal  expansion  tests,  and  simple  uniaxial  and  biaxial 
loading  configurations.  The  element  performed  properly,  passing 
all  tests. 

5.0  USER  INSTRUCT  [ONS  AND  SAMPLE  PROBLEMS 

The  following  sections  present  the  required  user  input  and 
detailed  instructions  for  using  the  developed  crack  element.  The 
first  problem  uses  a  fairly  crude  mesh  to  calculate  the  Kj  stress 
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Figure  9  Boundary  conditions  for  edge  crack  and  central  crack  specimens 


intensity  factor  for  an  edge  crack.  Another  more  involved  prob¬ 
lem  which  demonstrates  calculation  of  both  Kj  and  Kjj  is  also 
presented . 

5.1  REQUIRED  USER  INPUT  AND  INTERPRETATION  OF  OUTPUT 

The  crack  element  was  Implemented  using  the  CDUM1  user 
element.  The  formats  for  the  ADUM1,  CDUM1,  and  PDUM1  cards  are 
presented  in  Figures  11  through  13.  The  ADUM1  card  is  used  by 
NASTRAN  so  that  it  interprets  the  associated  PDUM1  cards  and 
CDUM1  cards  properly.  This  card  should  always  be  used  as  shown 
in  Figure  11  and  only  1  card  per  NASTRAN  bulk  data  deck  is 
required . 

The  PDUM1  cards  and  CDUM1  card  fields  are  described  in 
Figures  12  and  13.  Also  presented  in  Figures  12  and  13  are 
allowable  ranges  for  the  various  options. 

The  CDUM1  user  output  consists  of  nine  headings  labeled  SI 
through  S9.  Table  3  describes  the  various  output  quantities  for 
the  crack  element.  IKI  of  Table  3  is  input  on  the  element's 
PDUM1  card. 

Note  that  currently,  the  BANDIT  *  -1  option  must  be  used  on 
the  NASTRAN  card. 

5.2  SAMPLE  PROBLEM  FOR  CALCULATION  OF  Kj 

Figure  14  presents  a  crude  model  for  calculating  the  stress 
intensity  factors  for  the  edge  crack  with  uniform  stress.  This 
geometry  was  previously  depicted  in  Figure  7b,  with  the  associ¬ 
ated  errors  shown  in  Table  1  under  the  37  Grid  Mesh  column.  The 
bulk  data  input  is  shown  in  Figure  15. 

Elements  1  through  *1  are  singular  crack  elements  while 
Elements  5,  6,  and  9  are  non-singular.  Note  the  grid  numbering 
sequence  for  Elements  1  and  4.  Since  the  stress  intensity  factors 
are  calculated  along  rays  closest  to  the  element's  local  x-axis. 


BULK  DATA  DECK 


Input  Data  Card  ADUM1 


Dummy  Element  Attributes 


Description:  Defines  attributes  of  the  dummy  element  CDUM1 


Format  and  Example: 


ADUM1 


ADUM1 


Number  of  grid  points  connected  by  DUM1  dummy  element 
(Integer  =  8) 


Number  of  additional  entries  on  CDUM1  connection  card 
(Integer  =  1) 


Number  of  additional  entries  on  PDUM1  property  card 
(Integer  =  6) 


Number  of  displacement  components  at  each  grid  point 
used  in  generation  of  differential  stiffness  matrix 
(Integer  =  3) 


Figure  11  ADUM1  bulk  data  card  for  singular  or 
tion-o  i  ni'ular  structural  element- 


BULK  DATA  DECK 
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Input  Data  Card  CDUM1  Dummy  Element  Connection 

Description:  Defines  a  sinqular  or  non-sinqular,  two-dimensional  structural 
element.  Must  be  used  in  conjunction  with  the  Fortran  code 
supplied  with  this  report. 


Format  and  Example: 


1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

CDUM1 

EID 

PID 

Gl 

G2 

G3 

G4 

G5 

G6 

abc 

CDUM1 

114 

108 

2 

5 

6 

8 

7 

11 

ABC 

rbc 

G7 

G8 

X 

X 

X 

X 

X 

X 

+BC 

12 

14 

Field  Contents 


EID 


Element  identification  number  (Integer  >  0) 


PID  Identification  number  of  a  PDUMI  property  card 

(Integer  >  0) 

G1...G8  Grid  point  identification  numbers  of  connection  Doints 

( Interqer  0,  G1  f  G2  ...  f  G8) 


Remarks  1.  All  grid  points  must  be  unique  and  all  eiqht  qrids  must  be 

present.  Dummy  grid  points  must  be  introduced  into  the  model 
if  the  element  is  to  have  less  than  eight  grid  points. 

2.  To  form  a  trianqle,  the  x,y,z  coordinates  of  G4  must  equal 

the  x ,y,z  coordinates  of  Gl.  G4  is  then  SPC'd  in  123456  and  is 
considered  a  dummy  qrid  point  not  used  in  the  analysis. 

3.  To  eliminate  any  mid-side  qrid,  the  x,y,z  coordinates  of  the 
mid-side  qrid  must  equal  the  x,y,z  coordinates  of  one  of  the 
corner  grids.  The  eliminated  mid-side  qrid  is  then  SPC'd  in 
123456  and  is  considered  a  dummy  qrid  point  not  used  in  the 
analysis. 


--  continued 


Figure  1?  ODUM!  bulk  data  card  for  singular  or  non-alnrular 
rd.r-ud  iiral  element. 
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CDUM1  (continued) 


Remarks 


4.  Dummy  qrids  may  be  shared  by  adjacent  elements  in  order  to 
keep  down  the  total  number  of  required  grid  points. 

5.  The  ordering  convention  for  the  element's  grid  points  G1 
through  G8  are  shown  below. 


G4 


G8 1 


G7 
— - •- 


Figure  12  Concluded. 


BULK  DATA  DECK 


Input  Data  Card  PDUM1  Dummy  Element  Property 

Description:  Defines  the  properties  and  stress  evaluation  techniques  to  be 
Input  used  with  the  CDUM1  singular  or  non-singular  element.  Must  be 

used  in  conjunction  with  the  Fortran  code  supplied  with  this 
report. 


Format  and  Example: 


1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

PDUM1 

PID 

MID 

T 

GAMMA 

IPLANE 

NIPR 

NIPS 

IKI 

3DUM1 

108 

2 

0.10 

0.50 

1 

5 

4 

1 

Field  Contents 

PID  Property  identification  number  (Integer  >  0) 

MID  Material  identification  number  (Integer  >  0) 

T  Element  thickness  (Real  >  0) 

GAMMA  Exponent  used  in  displacement  field  (Real,  0.50  for 

singular  element,  2.0  for  non-singular  element) 

IPLANE  Plane  strain  or  plane  stress  option,  use  0  for  plane 

strain,  1  for  plane  stress  (Integer  0  or  1) 

NIPR  Number  of  integration  points  in  r  direction. 

The  r  direction  is  the  radial  direction  for  the 
singular  element  (0  <  Integer  <  5) 

NIPS  Numer  of  integration  points  in  s  direction 

(0  <  Integer  <  5) 

IKI  Stress  and  stress  intensity  factor  calculation  option. 

Use  I KI =0  to  calculate  K.  and  K.j  based  upon  stresses. 

Use  IKI=1  to  calculate  k|  and  K;!  based  upon  displacements. 
For  I K I =2  no  stress  intensity  factor  calculations  are 
performed.  (Integer  0,1  or  2) 


Figure  13  PDUM1  bulk  data  card  for  singular  or  non-sin. -,ul 
structural  element. 


ar 


30 


TABLE  3 


INTERPRETATION  OF  CDUM1  USER  ELEMENT  STRESS  OUTPUT 


IKI* 

SI 

S2 

S3 

S4 

S5 

S6 

S7 

S8 

S9 

0 

X 

y 

°x 

°y 

Txy 

KI 

KII 

0 

0 

1 

X 

y 

ax 

ay 

Txy 

KI 

KII 

0 

0 

2 

X 

y 

°x 

ay 

Txy 

ex 

ey 

Y 

xy 

0 

*  F 

'or  IKI 

=  0,  Kj 

and  Kjj 

based 

on  stresses 

(Equation  8) 

For  IKI  =  1,  Kj  and  K^  based  on  displacements  (Equation  9) 
For  IKI  =  2,  Kj  and  Kn  are  not  calculated 

NOTES:  x  and  y  are  the  element  coordinates  where  stresses 
and  strains  are  reported. 
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Figure  14  Sample  problem  for  calculation  of  K^. 
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NASTRAN  BANDIT* -1 

ID  KI  TEST 

APP  DISPLACEMENT 

DIAO  8 

TIME  15 

SOL  1 

CEND 

TITLE  -KITEST.NID  -  TEST  OP  CDUM1  ELEMENT,  SIDE  CRACKED  GEOMETRY,  6/10/85 
SUBTITLE*  EDOE  NOTCH,  37  NODES,  A/W-l/3,  SINGULAR  ELE,  5X4  OAUSS 

SPC-1 
LOAD-1 
DISP-ALL 
STRESS-ALL 
BEOIN  BULK 

t  1..  2..  3..  A..  5..  6..  7..  8.,  9..  10 

* 

$  PARAMETERS 


GRDPNT  0 
COUPMASS  1 


PQDMEM1 


-  1  SING  ELEM  W.  KI  AND  KII  BASED  ON  DISPLACEMENTS 
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Element  1  is  numbered  as  shown.  Since  the  Element  1  x-axis  edge 
forms  the  cracks  surface,  the  PDUM1  card  for  Element  1  should 
request  that  displacements  be  used  to  calculate  stress  intensity 
factors.  This  is  consistent  with  the  equations  presented  in 
Section  2.3.  Element  4  is  numbered  as  shown,  and  its  stress 
Intensity  factors  are  based  on  stresses.  Theoretically,  the 
stress  Intensity  factors  predicted  for  Element  1  and  Element  4 
should  be  the  same.  It  has  been  found  that,  usually,  the 
displacement-based  stress  intensity  factors  are  more  accurate. 

For  Elements  2  and  3,  no  stress  intensity  factor  calculations  are 
performed  since  they  would  not  be  meaningful.  The  numbering  for 
Elements  2  and  3  is  not  as  critical,  although  the  first  grid  must 
be  located  at  the  crack  tip  for  all  singular  elements.  Since 
element  stresses  are  in  terms  of  the  element's  coordinate  system, 
it  is  suggested  that  they  be  aligned  with  the  basic  system 
whenever  possible. 

Grids  50  and  51  are  dummy  grids  which  must  be  present.  Note 
that  they  may  be  shared  by  any  adjacent  elements.  In  order  to 
form  a  triangular  element,  the  x,y,z  coordinates  of  element 
grid  4  (G4)  must  be  the  same  as  those  of  element  grid  1  (Gl). 
Mid-side  grids  are  eliminated  by  making  their  coordinates  equal 
to  that  of  any  of  the  corner  grids.  Note  also  that  Elements  5 
and  6  share  the  dummy  grids  52,  53,  and  54,  while  Element  9  uses 
dummy  grids  55,  56,  and  57. 

The  theoretical  solution  for  this  problem  is  given  in 
Reference  2,  page  2.10  as 

Kj  =  o  P(|) 

=  1000  psl  ATI  Tin  *  1.786 

=  3165  lb-in-^/^ 


The  reported  values  from  the  NASTRAN  run  are: 

K  =  2965  lb-in-3/2  (Element  1) 

COD  =  2  x  4.397  x  lo-1*  in.  (grid  1) 

These  calculated  errors  are  slightly  different  than  those  reported 
in  Table  1.  Because  the  errors  reported  in  Table  1  were  calcu¬ 
lated  using  double  precision  arithmetic,  whereas  the  NASTRAN 
results  shown  above  use  single  precision  arithmetic  for  stress 
recovery . 

5.3  SAMPLE  PROBLEM  FOR  CALCULATION  OF  Kj  AND  Kn 

Figure  16  presents  the  fairly  detailed  model  used  to  calcu¬ 
late  Kj  and  Kjj  for  a  central  crack  in  shear.  This  geometry  was 
previously  depicted  in  Figure  7(3,  while  the  loading  and  boundary 
conditions  are  shown  in  Figure  10.  Figure  17  presents  the  bulk 
data  input  for  this  model.  Several  additional  subcases  are  also 
analyzed  in  the  bulk  data  of  Figure  17*  The  SPC=1  set  applies  to 
the  shear  problem  whereas  the  SPC=2  set  is  for  the  tension  prob¬ 
lem.  Other  input  data  are  documented  in  the  deck.  Figure  16 
also  shows  a  blown-up  view  of  the  left  side  of  the  crack  geometry 
so  that  it  may  be  easily  compared  to  the  bulk  data  cards. 

Grids  325  and  326  are  the  dummy  grids  for  crack  Elements  1001 
through  1008.  The  elements  used  to  model  the  right  side  of  the 
crack  are  Elements  2001  through  2008.  Other  CDUM1  Elements  used 
in  the  model  are  transition  elements. 

The  theoretical  stress  intensity  factors,  Kjj*  for  the  shear 
loading  is  given  in  Reference  2,  page  10.1,  as 
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Kjj  =  1000  psi  / n ( 1 ) In  *  1.04  =  1843  lb-in 


while  the  crack  opening  displacement  (COD)  is  not  reported, 


The  corresponding  Kjj  as  reported  by  NASTRAN  is: 

KII  =  1924  lb-in”3/2  (Subcase  4,  Element  1004) 

Again,  a  slight  discrepancy  exists  between  these  results  and 
those  reported  in  Table  2  due  to  the  single  versus  double 
precision  method  of  stress  calculations  previously  described, 


6.0  SUMMARY  AND  CONCLUSIONS 

The  theory,  implementation,  and  user  instructions  for  a 
linear  elastic,  two-dimensional  COSMIC/NASTRAN  crack  tip  element 
have  been  presented.  The  element  was  incorporated  using  the 
dummy  element  capability  of  NASTRAN.  Subroutines  for  calculating 
the  stiffness  matrix,  mass  matrix,  thermal  load  vector,  stresses, 
and  stress  intensity  factors  have  been  developed.  Instructions 
for  using  the  element  along  with  sample  problems  have  been 
presented . 

The  stress  intensities  and  crack  opening  displacements  ob¬ 
tained  using  the  element  have  been  compared  to  several  theoretical 
solutions.  Both  COD  and  stress  intensity  factors  appear  to  be 
accurately  represented  even  for  relatively  coarse  meshes.  The 
accuracies  obtained  are  well  within  the  accuracies  required  by 
typical  engineering  calculations,  since  the  scatter  alone,  in  the 

Kj  values  from  a  typical  test  may  be  10%. 

Finally,  due  to  the  generality  of  the  element,  extensions  to 
include  anisotropic  materials  and  plasticity  should  be  easily 
accommodated.  The  work  performed  here  provides  a  solid  theoreti¬ 
cal  and  implementational  basis  for  developing  an  enhanced  version 
of  this  element. 
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APPENDIX 


MAGNETIC  TAPE  FORMAT  AND  PROCEDURE  FOR 
INSTALLATION  OF  DUMMY  ELEMENT  CODE  ON  VAX  COMPUTER 

The  dummy  element  FORTRAN  source  files  and  sample  problems 
are  provided  on  magnetic  tape.  The  tape  is  written  in  VAX 
FILES-11  format  and  has  the  following  characteristics. 

1600  BPI 
9  TRACK 
ASCII 

To  install  the  code,  all  files  should  first  be  copied  from 
the  tape  to  either  a  user's  disk  space  or  to  the  system  NASTRAN 
account  disk  space. 

The  following  files,  in  order,  will  be  found  on  the  tape. 

Dummy  element  stiffness  and  mass  matrix  routines 
KDUM1.MIS 

Dummy  element  thermal  loads  routines 
DUM1.MIS 
EDTL.MIS 

Dummy  element  stress  routines 
SDUM11.MIS 
SDUM12.MIS 

Bulk  data  decks  for  sample  problems 
KITEST.NID 
K12TEST.NID 


To  install  the  subroutines,  they  must  be  compiled  and  linked 
to  the  existing  NASTRAN  executable.  Merge  all  the  FORTRAN  files 
(.MIS  extensions)  into  one  combined  file  called  DUM1.T0T.  The 
user  must  have  available  the  two  NASTRAN  libraries,  NASPRILIB.OLB 
and  NASSECLIB.OLB.  The  following  command  procedure  should  be 
stored  in  a  file  called  NC.COM  and  then  executed  by  typing  in  @NC 
DUM1.T0T. 

$SET  NOON 

$FORTRAN/NOF77  ' PI ' /OPTIMIZE/NOLIST/OB JECT=TEMPLIBR. OBJ 

$KLIB : ="NASPRILIB . OLB" 

$LEN  =  'F$LENGTH(P1 ) 

$IF  (  ,F$L0CATE(".VSS,',P1)  .NE.  LEN)  -  THEN 
KLIB: “"NASSECLIB.OLB" 

$LIBRARY/REPLACE/LOG  'KLIB'  TEMPLIBR.OBJ 

$DELETE  TEMPLIBR.OBJ;* 

$SET  ON 

The  above  procedure  compiles  the  source  and  inserts  it  into 
the  old  libraries.  This  procedure  (NC.COM)  and  those  that  follow 
(LINK05.COM,  etc.)  are  provided  by  COSMIC  on  the  tapes  containing 
the  NASTRAN  program.  Refer  to  the  "Supplemental  Documentation 
for  VAX  NASTRAN"  provided  by  COSMIC  for  further  elaboration. 

The  following  procedure  links  the  object  files  created  above 
and  creates  the  executable  NAST08.EXE,  The  following  should  be 
stored  on  a  file  called  LINK08.C0M  and  then  executed  by  typing  in 
6LINK08 . 

$LINK/NOMAP/EXE=NAST08 .EXE  NASPRILIB . OLB/INCLUDE* ( NAST08 ) , - 
NASSECLIB . OLB/ INCLUDE* ( VAX08 , vax08e , 1FTE2 ) , - 
NASPRILIB . OLB/LIB/INCLUDE*- 
(GPTABD,XSFABD,SEMDBD,TABFBD,CN36BD) 

$EXIT 


The  following  procedure  creates  the  executable  NAST05.EXE. 
The  following  should  be  stored  on  a  file  called  LINK05.C0M  and 
then  executed  by  typing  in  6LINK05. 

$LINK/NOMAP/EXE=NAST05 .EXE  NASPRILIB.OLB/INCLUDE«(NAST05) , 
NASSECLIB . 0LB/INCLUDE= ( VAX05 , I PTE 2 ) , - 
NASPRILIB.OLB/LIB/INCLUDE=- 
(GPTABD , XSPABD , SEMDBD, CN3  6BD ) 

$EXIT 

Finally,  the  following  procedure  creates  the  executable 
NAST13.EXE.  The  following  should  be  stored  on  a  file  called 
LINK13.COM  and  then  executed  by  typing  in  6LINK13. 

$LINK/N0MAP/EXE=NAST13.EXE  NASPRILIB.0LB/INCLUDE-(NAST13), 
NASSECLIB . OLB/INCLUDE* ( VAX1 3 , IFTE2 , SETC ) , - 
NASPRILIB.OLB/LIB/INCLUDE=- 

( GPT ABD , PLA4BD , XSPABD , SDR2BD , SEMDBD ,  -  CN3 6BD ) 

$EXIT 

After  executing  the  above  procedures,  three  new  files, 
NAST08.EXE,  NAST05.EXE,  and  NAST13.EXE  will  be  on  the  account. 
These  are  then  to  be  used  in  place  of  the  old  NAST08.EXE, 
NAST05.EXE,  and  NAST13.EXE  files  which  currently  reside  on  the 
NASTRAN  disk.  If  desired,  the  old  NAST08.EXE,  NAST05.EXE,  and 
NAST13.EXE  files  may  be  deleted  from  the  NASTRAN  disk  at  this 
time. 

The  Bulk  data  decks  KITEST.NID  and  K12TEST.NID  for  the 
sample  problems  are  used  in  exactly  the  same  manner  as  any  other 
bulk  data  deck. 
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